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Abstract 

New hybrid Molecular Dynamics-Monte Carlo methods are proposed to increase 
the efficiency of constant-pressure simulations. Two variations of the isobaric 
Molecular Dynamics component of the algorithms are considered. In the first, we 
use the extended-ensemble method of Andersen [H. C. Andersen J. Chem. Phys. 
72,2384 (1980)]. In the second, we arrive at a new constant-pressure Monte Carlo 
technique based on the reversible generalization of the weak— coupling barostat [H. 
J. C. Berendsen et al J. Chem. Phys. 81, 3684(1984)]. This latter technique 
turns out to be highly effective in equilibrating and maintaining a target pressure. 
It is superior to the extended-ensemble method, which in turn is superior to 
simple volume-rescaling algorithms. The efficiency of the proposed methods is 
demonstrated by studying two systems. The first is a simple Lennard- Jones fluid. 
The second is a mixture of polyethylene chains of 200 monomers. 

1 Introduction 

Molecular simulations at constant pressure are attractive because they mimic the con- 
ditions often encountered in laboratory experiments. This applies to both Molecular 
Dynamics and Monte Carlo calculations. In Monte Carlo simulations, trial configu- 
rations are generated in a random manner, and accepted according to well-defined 
probability criteria. The so-called "hybrid" schemes constitute one of many ways of 
proposing such trial configurations. In a hybrid Monte Carlo trial move, a short molec- 
ular dynamics trajectory is used to generate a new configuration the move is 
subsequently accepted or rejected according to probability criteria which take into 
account the way in which the trial configuration was generated (since molecules are 
moved in the direction of the forces acting on them, the resulting bias must be elim- 
inated). Molecular Dynamics becomes exact in the limit of vanishing time— step; the 
acceptance rate of trial configurations produced by a short MD run with small time- 
steps can therefore be expected to be high. The use of longer time-steps does not 
deteriorate the accuracy of a hybrid scheme, but it decreases the acceptance rate. 
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The Monte-Carlo acceptance/rejection procedure ensures that the correct statistical- 
mechanical ensemble is sampled, and the results remain exact, regardless of time-step. 
In contrast, a pure molecular dynamics simulation using large time-steps would lead 
to a fast algorithm but questionable results. 

Many algorithms for isobaric-isothermal NPT molecular dynamics have been de- 
veloped over the past decades {j-g. In contrast, NPT Monte Carlo simulations still 
rely on the use of simple, random volume changes accompanied by uniform rescaling 
of the molecules' coordinates. A few exceptions have been reported in the literature, 
but they lack the ease of implementation of simple random volume moves [Q. Un- 
fortunately, volume relaxation is notoriously slow in simple, random volume-rescaling 
Monte Carlo algorithms. Furthermore, their efficiency deteriorates as the size of a 
system increases. 

Our experience indicates that volume moves are actually one of the limiting steps 
towards successful Monte Carlo simulations of complex fluids. It has, for example, 
become clear in studies of solubilities of small molecules in polymers that the compu- 
tational bottleneck is indeed the equilibration of the density || . Interestingly, litera- 
ture results suggest that constant-pressure molecular dynamics algorithms can provide 
relatively fast volume relaxation. It is therefore of interest to explore whether a hybrid 
NPT method is able to circumvent some of the problems associated with conventional 
NPT Monte Carlo techniques. This paper considers two hybrid NPT approaches. We 
find both of these to be superior to conventional methods for constant-pressure Monte 
Carlo simulations. 



2 Constant Pressure Monte Carlo 

Before presenting the two hybrid methods to be discussed here, we shortly review the 
fundamentals of Monte Carlo simulations in the isothermal— isobaric ensemble. The 
natural variables are particle number N, pressure P, and temperature T; the statistical 
mechanical potential for this set of variables is the Gibbs free energy G(N, P,T). As 
N, P and T are fixed, their conjugate variables must fluctuate. These are the chemical 
potential /i (which, for a pure system, is equal to the Gibbs free energy density G/N), 
the system volume V and the entropy S. 

The relevant probability distribution function for an NPT ensemble is given by 

f NPT = exp[-(3U({r},0}) - pP*V({r})\- (1) 

Here P* is the system pressure, {p\ is the set of all momenta, {r} the set of all 
positions, U the internal energy of the system and V its volume; (3 = l/(fcsT), where 
ks is Boltzmann's constant and T is the temperature. 

This distribution function is to be sampled by a simulation in an NPT ensemble; 
the following acceptance criteria for trial moves generate configurations distributed 
according to Equation (|j): 

Pacc = min{l,exp[-/3A[/({r, {p}) - (3P*AV({r})]}, (2) 
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where A denotes the difference of a property between the new (proposed) state and the 
old (original) state, and where min(a;, y) denotes the minimum of its two arguments. 
The condition of detailed balance, namely 

K(a -> b) = K{b -> a) 

K(a ->b) = p(a)q(a -> b)p acc (a -> b), (3) 

facilitates the construction of a correct algorithm. In Equation (||), K(a —> b) is 
the probability of going from state a to state b; this probability can be written as 
the product of p(a), the probability being in state a, q(a — > b), the probability of 
proposing a transition from a to b, and p acc (a — > &), the probability of accepting such 
a transition. Accepting or rejecting trial moves according to 

*w(a -» 6) = mm 1, i£-L± -( (4) 

V P(a)q(a^b)J 

ensures that configurations are distributed according to the desired probability distri- 
bution p appearing in Equation (^). 



3 The extended ensemble algorithm 

In order for a Molecular Dynamics algorithm to be suitable for direct application within 
a hybrid Monte Carlo formalism, two conditions must be fulfilled. The algorithm 
should be reversible in time and symplectic, i.e. preserving phase— space volume. If 
these requirements are met, Equation (^) is fulfilled and the correct Markov process is 
produced. Given these constraints, the Verlet integration scheme has often been used 
in this context |J. It is reversible, it preserves phase— space volume, and it is correct 
to order 0(At 3 ). 

The extended-ensemble technique of Andersen || (in conjunction with the Verlet 
integrator) satisfies the above prerequisites. Its application in the context of a hybrid 
move is therefore promising. The relevant equations of motion are given by (for a 
derivation see the original paper of Andersen Q ) : 

d t n = ^ + ld t lnV 

N 

i=l i<j=l J 

where M is a fictitious mass of a "piston" associated with the volume move (the piston 
connects the system to an infinitely large pressure reservoir). Parameter a describes 
the fictitious potential of the volume. These equations can be easily implemented, 
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and the resulting molecular dynamics trajectory obeys a constant NPH constraint. 
Anderson's method can be generalized to constant NPT, but this is not of interest here 
as maintaining a constant temperature is left to the Monte Carlo acceptance criteria. 

In most Monte Carlo simulations one deals with particle positions only. Therefore, 
at the beginning of a hybrid step, velocities (or momenta) must be assigned to each 
particle. These are generated at random from a Maxwell-Boltzmann distribution. In 
addition, a "volume" velocity vy is also generated. 

The acceptance criteria for a volume hybrid move are derived from Equations (^|) 
and (^). Since the integration of the equations of motion itself is deterministic, after 
having assigned velocities to the molecules, the probability of reaching a final state 
from some original state is unity. Still, the probability of assigning velocities to the 
particles must be considered, thereby leading to a transition probabilities of the form: 

?(o-&) = Etlv{vf) 

q(b^a) = 4S(«y b) ) (6) 

where p(vy) is the probability of generating a volume velocity according to a specific 
distribution, Eki n is the kinetic energy, and superscripts a and b serve to denote the 
original and new (trial) states, respectively. The acceptance criteria for hybrid NPT 
moves are obtained by substituting Equations (Q) into Equation (^). 



4 The reversible Weak- Coupling algorithm 

One of the most efficient algorithms for constant temperature and/or pressure Molecu- 
lar Dynamics is the weak- coupling scheme proposed by Berendsen et al. ||. Although 
it is widely used, unlike Andersen's algorithm, it is not time-reversible due to the 
first-order nature of the equations of motion. This precludes its direct use in hybrid 
Monte— Carlo moves. The algorithm proposed in what follows eliminates this short- 
coming. 

Recently, it has been shown that the Weak-Coupling thermostat produces a correct 
ensemble, in agreement with other thermostat methods for any value of the coupling 
parameter r |u|. We briefly describe the Weak-Coupling algorithm, focusing only on 
the constant pressure case. In any Molecular Dynamics algorithm the integrator moves 
the velocities and positions in a time-step according to the equations of motion: 

n(t) ^ n(t + At) (7) 
Vi(t) -> Vi(t + At) (8) 

After this update, the weak-coupling algorithm calculates the new pressure P of the 
system based on the positions and velocities. This pressure is compared to the desired, 
target pressure P* . The volume of the simulation box is now scaled "towards" the 
target pressure, i. e. 

V(t + At) = r-'Vit) ■ [ P{t) p /* )At (9) 
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where t p is a characteristic time for this relaxation process (see below), and At is 
usually (although not necessarily) the same time-step used to integrate the equations 
of motion. The positions of all particles must be adapted to the new volume, i.e. 



The smaller the relaxation r p , the more closely the instantaneous pressure is tied 
to the target, and the stronger the disturbance of the actual dynamics by individual 
rescaling operations. In molecular dynamics, r p must be large enough to produce 
meaningful data, as the fluctuations of pressure produced by this algorithm are not 
entirely correct [jnj. This does not pose a problem for Monte Carlo, as the molecular 
dynamics moves are only used to propose trial configurations. 

As pointed out earlier, the volume rescaling of the Weak-Coupling algorithm is 
not time-reversible. This violates the detailed balance rate equation (^|). The solution 
to this problem is to allow the algorithm to run explicitly backward and forward in 
time; at the beginning of a hybrid move, a random number £ is drawn between and 
1. If £ < 0.5, the time direction in the simulation runs forward; otherwise we run 
backwards. Running the simulation in reverse is equivalent to setting the relaxation 
time t p to a negative value, as the rest of the equations of motion is symmetric under 
time-reversal. The acceptance criteria for this constant pressure hybrid move is similar 
to that for the Andersen case. 

The weak-coupling method as outlined so far results in rather slow fluctuations 
around the desired pressure. Furthermore, (if backward moves are momentarily dis- 
regarded) the relaxation becomes exponentially slow close to the target pressure. 
This limits the efficiency of the simulation. To improve efficiency, we generalize 
the weak-coupling hybrid scheme to a target pressure selected randomly within a 
bounded domain. The pressure P* used as local target pressure for the simula- 
tion in one hybrid step is now chosen with uniform probability from the interval 
[(1 — -P ra nge)-f^ ext \ (1 + -Prange)-P^ cxt ^] , whereas the pressure in the acceptance criteria 
remains untouched (at p( cxt )). Values for P ra ngo range from to 0.1. Since the local 
target pressure is uniformly chosen from the neighborhood of the global target pressure, 
detailed balance is still obeyed. 

5 Efficiency Considerations 

The main parameter which can be tuned and optimized in the reversible weak-coupling 
(RWC) simulation is the relaxation time t p . The smaller r p , the closer the system 
follows the desired target pressure, but the lower the acceptance rate is. 

We analyze the density autocorrelation function to provide a measure of efficiency 
for the proposed hybrid schemes. Figure 0a) shows that, for a simple Lennard Jones 
fluid, increasing the time-step leads to smaller correlation times and faster equilibration 
of the density. However, if the ratio of correlation time to time step becomes too small 
the simulation becomes again inefficient as one would need extremely large neighbor 
lists (or extremely frequent updates) to balance the rapid motion. The figure shows 
how important a good choice of the parameters can be. An unfortunate choice results 
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Monte Carlo Steps 

Figure 1: Density Autocorrelation function for a Lennard- Jones fluid (T = 1.4, p = 
0.21732, — > p = 0.45) with different values of time-step and correlation time (part a). 
Prange = 0.1 the values for the timestep and the correlation time are given in the graph. 
Part b) shows the comparison to simple rescaling and the hybrid method based on the 
Andersen barostat. The box mass M in the Andersen case was set to M = 1.0. 

in only marginal improvements relative to conventional methods. 

The ratio between the time-step and the correlation time is important. The cor- 
relation function for At = 0.001, r = 0.25 is the same as for At = 0.002, t = 0.5 (not 
shown). We found that a value of this ratio of about 250 provides good efficiencies. If 
we increase this ratio the acceptance rate drops; if we decrease it the efficiency of the 
single move deteriorates. 

The conventional method for constant pressure Monte-Carlo simulations comprises 
simple volume-rescaling trial moves. In such moves a random change of volume is 
proposed according to 

Kcw = Kid + AV, (10) 

with Ay chosen uniformly from [—AV m&K , AV max \. 

Our results indicate that, even for a simple Lennard- Jones fluid, the hybrid tech- 
niques proposed here provide methods that appear to be superior to the simple volume- 
rescaling technique. For a good choice of parameters, the efficiency gains can amount 
to more than an order of magnitude. 

6 Test cases 

6.1 The Lennard- Jones Fluid 

As a test case, the Lennard- Jones fluid is revisited. We simulate a system of 200 
Lennard-Jones particles at a range of temperature and pressures. The cutoff is r = 2.5a 
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-5.27 



Table I: Comparison of density and internal energy per particle for the Lennard-Joncs 
system from simulation of this work using reversible weak coupling (1) and from the 
work of Johnson et al. (2). 



P KTl K-T2 

013 6.7 ±2 — 

0.40 2.4 ±0.5 2 

0.84 0.04 ±0.01 0.05 

Table II: Isothermal compressibilities for the Lennard- Jones fluid at T = 1.6 calculated 
by fluctuations of the volume (kti) an d by the equation of state by Rosenfeld (kt2)- 



and the usual analytic long-range corrections for energy and pressure are applied |12] . 
All values are given in the standard dimensionless units. The simulations used 1 million 
MC steps, each containing 5 constant-pressure molecular dynamics integration steps. 

We compare our results to those obtained from an equation of state jl3| which is 
known to provide accurate properties for the Lennard- Jones fluid. At a given tem- 
perature and pressure we compare density and internal energy. The results are shown 
in table ^ agreement between simulations and equation of state predictions is excel- 
lent, showing that the algorithms proposed here lead to correct property predictions. 
Note that we simulated the same systems (using identical starting conditions) using 
the three methods considered in this work (Andersen, Weak-Coupling, and Simple 
Volume- Rescaling). There were no noticeable differences in the results. 

We now examine the volume fluctuations produced by the different methods. These 
fluctuations are of interest because they provide a stringent test of the correctness and 
the efficiency of a simulation technique, and they are related to the compressibility of 
the fluid. We compare our results to the equation of state by Rosenfeld |l4[ , again at 
T = 1.6 (cf. Table ||). As expected, the compressibility corresponding to our various 
simulation techniques agrees with accepted values. 

KT= k B TV ■ 
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Figure 2: a) Density and b) Pressure as a function of simulation time for the reversible 
Weak-Coupling algorithm an the volume-rescaling algorithm for a system of 11 alkanc 
chains of length 200 and 26 ethylene molecules at P = 260kPa, T = 400X. For the 
pressure, a running average over 50 data pointsclarity. 



6.2 Polymer Melts 

In order to examine the performance of the proposed hybrid methods in the context 
of a complex fluid, we used them to equilibrate a polymer melt. 

Systems consisting of linear polyethylene chains (Cioo to C500) and a small frac- 
tion of ethylene were placed in a simulation box at a density above the experimental 
value. Three identical systems were left to equilibrate by the reversible weak-coupling 
hybrid methods and the volume-rescaling algorithm. The weak-coupling algorithm 
proved to be efficient and produced reliable, well-equilibrated densities for the pres- 
sure range under study (200 — 500kPa). All studies have been performed using the 
NERD force-field |,[l|. 

Figure shows a comparison between the convergence of the newly proposed weak- 
coupling method and simple volume rescaling. Both simulations started at a system 
density of 0.82 g/cm 3 , which is slightly too high. The equilibrium density is around 
0.75 g/cm 3 . It can be seen in the figure that the convergence of the weak-coupling 
algorithm to the correct, target pressure is relatively fast. In contrast, pressure (or 
volume) relaxation in t he conventional volume-rescaling simulation is slow. Even 
after 500, 000 steps, the running-average pressure in the conventional volume-rescaling 
simulation is far from the correct, target result (by about a factor of 3). 



7 Conclusions 

New hybrid schemes have been presented for NPT Monte Carlo simulations. We find 
that the proposed methods are more effective at relaxing the volume of simple and 
complex liquids than conventional Monte Carlo techniques. The increase of efficiency 
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becomes more apparent in fluids that exhibit long relaxation times (e.g. polymers). 

An additional advantage of the methodologies discussed in this work is that avail- 
able simulation codes for molecular dynamics calculations can be modified to perform 
exact Monte Carlo calculations with minimal changes. Implementation of the proposed 
techniques does not add a significant overhead, and the new moves not only change 
the volume, but also contribute significantly to the repositioning of particles, thereby 
improving sampling. 

It is our experience that equilibrating the density of a dense polymer melt using 
conventional NPT Monte Carlo moves is highly computationally demanding, and often 
not possible at all. In fact, we have shown in this work that, as a result of poor 
sampling, the compressibility of a melt of short polyethylene chains determined using 
simple rescaling techniques deviates considerably from the correct value. In contrast, 
the compressibility predicted using the weak-coupling algorithm is in agreement with 
experimental values. 
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